Contents

Wannierization

Electronic simulation codes like Quantum Espresso calculate the Bloch bands \(\varepsilon_{\mathbf kn}\) and wavefunctions \(\psi_{n\mathbf k}(\mathbf r)\) of a crystal lattice. These bands are gauge invariant while the wavefunctions allow for an ambiguity in the choice of complex phase. These gauges need to be consistent if we wish to transform these wavefunctions to a localized Wannier function [MMY+12]:

(3)\[ |\mathbf Rn\rangle=\frac{V}{(2\pi)^2}\int_{BZ}\text d\mathbf k\, \text e^{-i\mathbf k\cdot\mathbf R}|\psi_{n\mathbf k}\rangle. \]

Projection method

One way to remove the gauge inconsistency is to project trial orbitals \(|g_n\rangle\) on the basis of the Bloch wavefunctions:

(4)\[ |\phi_{n\mathbf k}\rangle=\underbrace{\sum_m|\psi_{m\mathbf k}\rangle\langle\psi_{m\mathbf k}}_{P_\mathbf k}|g_{n}\rangle. \]

The projector \(P_\mathbf k\) is gauge invariant so the \(|\phi_{n\mathbf k}\rangle\) are as well. The \(|\phi_{n\mathbf k}\rangle\) are however not necessarily orthogonal. In general we have:

(5)\[ \mathbf S_{mn,\mathbf k} = \langle\phi_{m\mathbf k}|\phi_{n\mathbf k}\rangle\neq\delta_{mn} \]

Marzari et al. [MMY+12] propose to orthogonalize the states using a Löwdin transform which finds a set of orthonormal basis functions \(|\tilde\psi_{n\mathbf k}\rangle\) that minimizes:

(6)\[ \sum_n \int\text d \mathbf x\,\left |\phi_{n\mathbf k}(\mathbf x)-\tilde\psi_{n\mathbf k}(\mathbf x)\right|^2, \]

which corresponds to taking:

(7)\[ |\tilde\psi_{n\mathbf k}\rangle=\sum_m|\phi_{m\mathbf k}\rangle \left(\mathbf S_\mathbf k^{-1/2}\right)_{mn}. \]

Example: 1 band

We illustrate the projection method using an example of a 1d isolated band with mock Bloch wavefunction \(\psi_\mathbf k(\mathbf r)=u_\mathbf k(\mathbf r)\text e^{i\mathbf k\cdot\mathbf r}\). This functions must obey the periodicity requirements:

  • \(\psi_{\mathbf k+\mathbf G}(\mathbf r)=\psi_{\mathbf k}(\mathbf r)\), with \(\mathbf G\) a reciprocal lattice vector \(\rightarrow u_{\mathbf k+\mathbf G}(\mathbf r)=\text e^{-i\mathbf G\cdot\mathbf r}u_{\mathbf k}(\mathbf r)\).

  • \(u_\mathbf k(\mathbf r+\mathbf T)=u_\mathbf k(\mathbf r)\), with \(\mathbf T\) a lattice vector.

We take lattice spacing \(a=1\) in 1d with mock Bloch wavefunction:

(8)\[ u_{k}(x)=\cos(\pi x) (\cos(6\pi x)-2\cos(k))\cdot \text e^{ixf(k)}, \]

../_images/dft_2_0.png

Fig. 12 Phase of \(u_k\) compared to linear dependance on \(k\).

with \(f(k+2\pi)=f(k)+2\pi\):

\[ f(k) = \pi\left(1-\cos^5\frac{\text{mod}(k, 2\pi)}2\right) + (k - \text{mod}(k, 2\pi)) \]

We plot \(\psi_{k}(x)\):

We note that the real part of \(\psi_k(x)\) is symmetric while the imaginary part is antisymmetric. Any symmetric \(g\) will then be purely real. For this function we choose a Gaussian:

(9)\[ g(x)=\text e^{-x^2}. \]

Any normalization of \(g\) and \(\psi\) will cancel in the construction of the orthonormal \(\tilde\psi\). We now do the inverse transform of the Bloch wavefunctions of ref{eq:wannier_transform} do get the following localized Wannier function:

We see that the Wannier function is localized around \(x=0\) with local maxima (minima) appearing around \(x\in \mathbb Z\) which become smaller for larger \(x\). We explain this behavior by taking some \(u(x)\) independent of \(k\). Now the Wannier function for a single band is quite simple:

../_images/dft_7_0.png

Fig. 13 The sinc function which forces the Wannier orbital do decay with a power law for a single band with \(u\) independant of \(k\).

(10)\[ |\mathbf Rn\rangle=\frac{V}{(2\pi)^2}u(\mathbf r)\int_{BZ}\text d\mathbf k\, \text e^{i\mathbf k\cdot(\mathbf r-\mathbf R)}\propto u(\mathbf r)\text{sinc}(\pi(\mathbf r-\mathbf R)), \]

which means that any wannierization of a single band boils down to removing the complex phase and using a sinc function as a weight to restrict the \(u\) to a unit cell as best as possible. This also means any decay of a Wannier orbital with \(u\) independent of \(k\) follows a power law.

Example: 2d atomic orbitals as Bloch wavefunction

To make matters a bit more interesting we’ll look at the hibridization of the s, p_\(x\) and p\(_y\) orbitals in a triangular lattice. As \(u\) component of the Bloch waves we take the three sp\(^2\) orbitals which we label with \(\gamma_i\) and as \(g_n\) we take the normal orbitals. Since we need to adhere to periodicity of \(u\) we take:

(11)\[ u_i(\mathbf r)=\sum_\mathbf R \gamma_i(\mathbf r+\mathbf R), \]

with \(\{\mathbf R\}=\left\{\left(0, 0\right), \left(\pm 1,0\right), \left(\pm 1/2,\pm\sqrt{3}/2\right), \left(0, \pm \sqrt3\right), \left(\pm 3/2, \pm\sqrt3/2\right)\right\}\), integration beyond NNN is exponentially surpressed.

../_images/dft_11_0.png

Fig. 14 Three periodic \(u_i(\mathbf r)\) in a triangular lattice based on sp\(^2\) orbitals. We aim to construct localized Wannier orbitals using the three standard atomic orbitals: s, p\(_x\) and p\(_y\)

Since calculation of integrals in 2d is much more computationally exhausting we move to a mesh representation of the functions. Any integration is just the sum of the mesh times the mesh tile size. We write ref{eq:nonorthogonal} out in matrix form:

\[\begin{split} \begin{pmatrix} |\phi_{1\mathbf k}\rangle\\ |\phi_{2\mathbf k}\rangle\\ |\phi_{3\mathbf k}\rangle\\ \end{pmatrix}=\underbrace{ \begin{pmatrix} \langle u_1|\text e^{-i\mathbf k\cdot\mathbf r}|g_{1}\rangle & \langle u_2|\text e^{-i\mathbf k\cdot\mathbf r}|g_{1}\rangle &\langle u_3|\text e^{-i\mathbf k\cdot\mathbf r}|g_{1}\rangle\\ \langle u_1|\text e^{-i\mathbf k\cdot\mathbf r}|g_{2}\rangle & \langle u_2|\text e^{-i\mathbf k\cdot\mathbf r}|g_{2}\rangle &\langle u_3|\text e^{-i\mathbf k\cdot\mathbf r}|g_{2}\rangle\\ \langle u_1|\text e^{-i\mathbf k\cdot\mathbf r}|g_{3}\rangle & \langle u_2|\text e^{-i\mathbf k\cdot\mathbf r}|g_{3}\rangle &\langle u_3|\text e^{-i\mathbf k\cdot\mathbf r}|g_{3}\rangle\\ \end{pmatrix}}_{A_\mathbf k} \begin{pmatrix} e^{i\mathbf k\cdot\mathbf r}| u_{1}\rangle\\ e^{i\mathbf k\cdot\mathbf r}| u_{2}\rangle\\ e^{i\mathbf k\cdot\mathbf r}| u_{3}\rangle\\ \end{pmatrix}, \end{split}\]

then \(S_\mathbf k=A_\mathbf k^\dagger A_\mathbf k\), which we diagonalize: \(S_\mathbf k=U_\mathbf k^\dagger S^d_\mathbf k U_\mathbf k\). Now we get: \(S_\mathbf k^{-1/2}=U_\mathbf k \left(S^d_\mathbf k\right)^{-1/2} U_\mathbf k^\dagger\) and \(S_\mathbf k^{-1/2}A_\mathbf k\) is the rotation matrix for the Bloch functions. We remind that any normalization of \(g_i\) cancels in this product.

../_images/dft_14_0.png

Fig. 15 The three Wannierfunctions which resemble the \(g_i\): \(\tilde\psi_1\) takes after \(p_x\), \(\tilde\psi_2\) after \(p_y\) and \(\tilde\psi_3\) after \(s\).

The three Wannier wavefunctions in The three Wannierfunctions which resemble the g_i: \tilde\psi_1 takes after p_x, \tilde\psi_2 after p_y and \tilde\psi_3 after s. resemble the atomic orbitals chosen for the \(g_i\). The \(p\) orbitals delocalize to neighboring sites.

Supercell Code